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ABSTRACT 

The details are given of a computer program which calculates the 
minimum time required for a space vehicle to meet its mission while in 
an optimum trajectory. Also, the program can maximize an option if the 
burning time is fixed. The vehicle motion, assumed to be in a plane, is 
described by two linear differential equations. A precoded subroutine, 
Runge-Kutta Fourth Order, is used to integrate the equations of motion. 
Numerical methods are employed to determine the deviation in constraining 
functions explicitly in terms of first and second order variations in 
thrust angle. An iteration method is employed which, when computed, 
provides the effect of first and second order variations in addition to 
the usual numerical solution to the differential equations. 
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DEFINITION OF SYMBOLS 




i 


* 


Symbol 

x 

y 

> 

X 

y 

X 


y 

V 

r 

e 


go 


Ti 

fc k 


A t k 


Definition 

Horizontal component of acceleration 
Vertical component of acceleration 
Horizontal component of velocity 
Vertical component of velocity 
Horizontal component of distance 
Vertical component of distance 
Velocity of vehicle 

Distance of vehicle from center of earth 
Angle of velocity vector from the vertical 
Angle of departure 

Gravity acceleration at initial time 

Distance from center of earth at initial time 

Horizontal drag component 

Vertical drag component 

Staging times 

Integration times 

Integration time intervals 




A £ 


k 


* 


fc k + tk+i 
2 




k-L 


v 


DEFINITION OF SYMBOLS 
(Continued) 

Definition 


fc k - T i 


Constant times force 


Constant times mass 


Constant times rate of change of mass 


Thrust acceleration 


Thrust angles 

Thrust angle at initial time 
Deviation of thrust angles 
Thrust angle at cutoff 


Orbital velocity 
Escape velocity 


Angular momentum 


Radius 90 from perigee 


vi 
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COMPUTER PROGRAM TRAJECTORY OPTIMIZATION 
BY EXPLICIT NUMERICAL METHODS 

SUMMARY 

The details are given of a computer program which calculates the 
minimum time required for a space vehicle to meet its mission while in 
an optimum trajectory. Also, the program can maximize an option if the 
burning time is fixed. The vehicle motion, assumed to be in a plane, is 
described by two linear differential equations. A precoded subroutine, 
Runge-Kutta Fourth Order, is used to integrate the equations of motion. 
Numerical methods are employed to determine the deviation in constraining 
functions explicitly in terms of first and second order variations in 
thrust angle. An iteration method is employed which, when computed, 
provides the effect of first and second order variations in addition to 

•A- 

the usual numerical solution to the differential equations. 


Vc 

See reference 1. 



INTRODUCTION 


The purpose of this report is to document a computer program 
"Trajectory Optimization by Explicit Numerical Methods," The program 
contains thirteen options with a maximum of eight additional constraints 
for each option and is applicable to space vehicles with one or more 
stages. There are three precoded subroutines that must be contained in 
the source deck." Three select switches, SS 2, SS 3, and SS 4 are 
used for additional printed data. SS 2 is turned on if the conditions 
of power cutoff are desired. If SS 3 is turned on, data throughout the 
trajectory is printed for all iterations. If it is off, only the last 
iteration is printed. Turn SS 4 on for the data at every integration 
interval during the coast period. No disks, drums, or special tapes 
are used. The program language is FORTRAN IV and the computer is 
Control Data Corporation's 3200. 

"Trajectory Optimization by Explicit Numerical Methods" was conceived 
and derived by Lyle R. Dickey, Technical and Scientific Staff, Aero- 
Astrodynamics Laboratory, George C. Marshall Space Flight Center. 


* See Appendices I, II, and III. 
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ANALYSIS AND EQUATIONS 




* 


* 


The results desired require the solution to the linearized 
equations which describe the motion of the vehicle. For analysis, the 
motion was assumed to be in a plane and described by the following 
differential equations : 


X = 4 Sin ^ + ’ 


y = — costf + y 

J m ; g 


where 


Xg= — g 


y 

y g = i g 


g = -g r 




e = 9 81 meters 
S ° second y ’ 


r G = 6370 KM , 


r = )f x^-f y 2 , 


v 




The analysis for — and are discussed on pages 5 and 6. 


* 
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At each staging time a b^, and are given for i=l,2,3, 

— >N. 

The numerical method, Runge-Kutta Fourth Order, is used for 
integrating the equations of motion. DELT, the integration step size, 
and are used to compute two tables of time, t^ and 



fc k = **-1+ DELT > 


At k = DELT 

for 

T i“ l k - DELT - 

If 

T.- t k <DELT , 

then 

" T i 


t k+l = T i ’ 


At k = T i- V 


At k+1 = °-°- 

For the case where 

T • — T . , < DELT 
i l-l 

then 

C k = T i-i > 


C k+1 = T i-1 » 


fc k+2 = T i-i + \ ( T i — T i_l^ » 


4 


fc k+3 " T i ’ 


fc k+4 T i ’ 


A t k = 0.0 , 


At k+i ■ | < T i - T i-i> > 


At k+2 t k+L ’ 


At k+3 " 0 - 0, 


The end points are required to be the same. 


fc i = T i • 


fc n T N 


Briefly stated, the t table is discontinuous at each staging time. 

An added input value of time, T2, is used for calculating the t table 
into irrational numbers. 


The thrust angles, are computed as follows: 


f k =i 0 + ( ^n " ^o> 


ti - 


k"ti 




where 


f Q = tan 


-1 


x 0 yo-yo x o 

x o*o + yoy 0 


‘See Deck Setup Instructions for choosing a proper value for T2 
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and x Q , y Q , x Q and y Q are initial conditions for integration. 

Thrust, ^ , is described by the following: 
f 

m k 

where T 

k 

Runge-Kutta requires the f angles to be computed once at t k , 
twice at t^+ I and once at t k+ ^ . Three point interpolation is 

used to compute ^ at t^+ At^ . 

t = t k + 2 At k 5 

c = ^“tk+l) ( t-t k) 

1 ^k-l-tk+l^tk-l-tk) 

C = *~k+j) ( tk_ i ) 

2 (tk-tk+iHtk-tk-i) 5 

C 3 = (t ~ t k^ (t " t: k-l) 

(t k+l _t k) ^k+l^k-l) 

then i, = C 1 ^ k _ 1 +C 2 tk+C 3 ^k +1 • 


b i +C i\ 


tk-Tj 

100 


For iteration there are six other tables computed at t. + |a t. 



4 _ t k -l+t k 

- 2 


^k £k-l > 


M = r] 

[^] • 


t- §o r o 

h 2 - ~ F* 


K, = h. 


K 0 = - 


go r o £ 


[- <i)1 


The following matrix is defined for iteration: 


«vW ‘ I+A< 4^ t k > 


for 


and 


k = 1 ,2,3, ,n, 


A(^) = 


0 

0 


K, K n 


1 

0 

0 

0 


The cutoff options and additional constraints are defined and 
stored in a single array, Z. 


Option 1. 


2 . 


Zi = Time = t k 
Z 2 = x = /f(x)dt 
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X 

C 4 = 

[ 0 

0 

1 

0 ] 


y 

C 5 = 

[ 0 

0 

0 

1 ] 

# 

r 

C 6 " 

[ - 
l r 

1 

r 

0 

0 ] 

1 

V 

C 7 " 

[ 0 

0 

X 

r 

! ] 


e 

C 8 = 

r xvcos0-xr 
L r^vsinQ 

yvcos0-xv 

r^vsinS 

xrcos0-xv 

rv^sinQ 

yrcosS-yv 

rv^sinfi* 


$ 

C 9 = 

1 

» — 1 

X 

~ P 

0 

0 ] 


Oy 

C 10 = 

r §o r o 2x 

1 3 

L r 

§o r o 2 y 

2x 

2 y ] 

* 

E v 

C U= 

r §o r o 2 x 

1“ r 3 

So r o 2 y 

pr- 

X 

y ] 


E 

c 12“ 

[-y 

• 

X 

y 

-X ] 


E. 

J 

c 13 = 

r x , 2y z 12 

+ p 

r §o r o 

y 2xZ i2 

r So r o 2 

2 y z i 2 

g r 2 
6 o 0 

2xZi2 

a r 2 J 

&o 0 


and 


-pm 


= Z 


pm 


"qm 


= Z 


qm 


♦ 


* 


where p and q are the numbers that specity power cutoff and 
coast cutoff options respectively for m=l,2 5 3 5 5 M 2 . 
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Kp m and Kp m are the values of power and coast cutoff limits for each 
option for m=l ,2,3, ,M 2 . 


T is defined when Z pm = K pm , 


at this time set — = 0.0 and integrate to 


T" , defined when Z ' = K“ 

m ? pm pm 


TLp m and TLp m are tolerances for the options. 


p* = C qm( T m ) 


m C qm (T*)x(T*) 


where 


*fl£> = 


y n 


y n 


The additional constraints are chosen such that 


P < r mJ> ' C Vj< T m> [ I+ 4 < T m> P m] 


for j =1,2,3, , N 


m 


then 


= 

m 


P < r n, £ > 

P <W 


p 'Vl ) 


m 
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p -c-oft 


(B) 


then 


Q - P M * 


or 


Q = P . 


If P is used, then q must be the number which specifies the 
time option, q = 1. 


where % = ~ (t) 


hj_ (t) = QU(T M “ ,t)H 1 (t) , 

h 2 (t) = QU(T M * ,t)H 2 (t) 

0 
0 

co s^/ 

— sini/ 


and K - |( * ) 


0 

0 

— sini^ 

— COS if 


F m (t) = 
m' ' 


C o St ^ T 
0.0 , T<t<T M " 
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/F(t)F'(t)dt = J 0 F(t)F' (t££ k , 

: o 

^ n 

- /F(t)h 1 (t)dt = (t)h x 

c o 

-l— 

B Q , 


= hi (t) + A'F(t), 



h* (t) = h 2 (t) +A'G(t). 


Zr u Of) 
Z r i 2 ( T 1 ) 


Zr lN x < T l> 

Z r 2i (Is ) 

^ r 22 (*^2 ) 

• • i 

Zr M 2 NM 2 < T M 2 ) 


then Az = K^— Z 


where = 


K, 


K 0 


%M 


K L5 given by the user, are the limits ot 
constraints . 


the additional 
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for 

user 


then 


I 6 (t k )| <E q , (C) 

I Az| < E1 l (d) 

=1,2,3, ,N^. E q and El^ are tolerances given by the 


If conditions C and D are not satisfied, then, 






5 


_ _ 

B = /F(t)F'(t)dt = 2 F(t k )F'(t k ^4 k , 
*b 

^ n 

£ „ -/nt)?(Odt - k i 1 r(c t )€(t k )4e k , 


AA 







= - e (t) +AA'F(t) . 


V' = i Jr + A , 
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return to integration with all the initial conditions and repeat 
the iteration until conditions C and D are satisfied. 



JOB NUMBER 560490 


DOCUMENTATION CHECKLIST 


Project: Have It 


1125 X 

Program Development Charts X 

Original Problem from Customer X 

Correspondence X 

Check Case X 

Documentation Checklist X 


Problem: 

^Abstract X 

''Technical Description X 

Equations X 

Definition of Terms X 

Special Options X 

Numerical Methods for Solution X 

Technical References X 

Related Projects X 


Programming: ■ 

Library Subroutines X_ 

Program Subroutines X_ 

Special Input Form X_ 

Special Input Tape Contents 

Special Output Tape Contents 

Plots Generated 

Method of Verification X 

*Block Diagram X 

Flow Charts X 

^Program Listing X 


Remarks 


None Used 
None Used 
None 


Deck Setup: 


^Computer Configuration X_ 

''Estimated Running Time X_ 

^Restart Procedure 

*Deck Sequence X_ 

*Required Input Data X 

^Restrictions 6c Limitations X 

^Diagnostics 

^Normal Output Quantity X 

^Definition of Output X 

^Sample Operator Instruction Card X 

^Sample Save Label 

*S ample Input Form X 

^Sample Output X 

* Sample Plots 


Decks : 

Symbolic 

Binary 

Data 


X_ 

X_ 

X 


No Capability 


None Printed 


Not Applicable 


None Used 


Location of 

Bldg. 4200, Rm. G24 
Bldg. 4200, Rm. G24 
Bldg. 4200, Rm. G24 


Prepared by: James W. Hilliard 

Programmer Date 


Approved by: Joseph Scollard 

Programming Mgr. Date 
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DECK SETUP INSTRUCTIONS 


All the input data must appear sequentially in the data deck as 
shown in Table 3. The user may not be concerned with T^ that must 
satisfy Equation A. T^, Ajg, B^, and Cjj must be added for all cases. 

T n <T l + DELT [175-2 (NfNO) ] , 

A N = A N-1 ’ 

% = ®N-1 » 

C N = C N-1 » 

Where N=1 + number of stages and 175 is the largest allowable number 
of storages for each single array variable. 

NO is the number of options and NR is the maximum number of 

constraints being used. NC=0 if Equation B is not used and will equal 
^ /\ 

q for all other cases. The value of q will be the number of the 
variable being maximized. 10 is the number of the cutoff option for 

which the weighting functions will be computed. 10=1,2,3, ,8. 

NRUN=0 for the last case in a stack. EO is the tolerance for €(t^.). 
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T2 is the second value of the t table desired to irrational ize all 


the values. 


and 


T2 = t 2 

T2 - T-l <DELT. 


The numbers for power cutoff, N6(I), coast cutoff options, N5(I), 
and the constraints, N2(I,J), can be obtained from the Z array. 

Refer to the Example Case and Data Package Card Arrangement, 

Table 4, for data instructions of all variables. 
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SAMPLE INPUT 


x o = 

236255.46 

y 0 = 

6456667.82 

*o = 

2800.5746 

y 0 = 

635.94317 

T i - 

184.99481 

T 2 = 

412.5503 

T 3 = 

523.87973 

T 4 " 

528.87973 

a i = 

10.787315 

a 2 = 

8.92405155 

a 3 - 

io- 20 

a = 

4 

1.96132994 

b i » 

1.2867664 

b 2 = 

0.69140715 

b 3 = 

1.0 

b 4 " 

0.36021199 

C 1 = 

-0.261282644 

C 2 = 

-0.212616813 

C 3 = 

0.0 

C 4 = 

0 . 046948359 

A t = 

5.0 

(1 

< 

5.0 
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SAMPLE INPUT 
(Continued) 


/\ 

q = 

7 


r o = 

6370000.0 


So = 

9.81 


e c = 

1.5707963 


*i = 

412.85503 

for time (t^.) 

= 

773.599673 

for time (t^) 

< - 

6370000.0 

for radius (r^.) 

K 2 = 

773.599673 

for time (t^) 

Pi = 

1 


P 2 = 

1 


§1 = 

6 


g 2 = 

1 


TL, = 

0.01 


TL„ = 

0.01 


TL“ = 

0.1 


TLj = 

0.01 


E o " 

10" 5 



5971137.6 

constraint for (y^) 

II 

H 

CM 

6552000.0 

constraint for (r^) 

K 

2; 2 

1.5707963 

constraint for (0^) 

EL = 
1A 

0.1 


E %- 

0.1 


E V 

io- 5 



This is a four stage vehicle, because there are four T points. 


N =1 + 4. 


There are two options for power cutoff and coast cutoff. 


NO = 2. 


The maximum number of additional constraints used is two. 


NR = 2. 


Equation E is used because q is given and NC = q if given, 


NC = 0. 


NC = 7 


The last option is used to compute the weighting functions 


10 = 2 


There is only one case. 


NRUN = 0 . 

E0 = 1.0 x 10 


T2 * 185.0. 


otherwise 


From the Z array the time option is 1 and the radius option is 6. 



N2 (1 , 1) 


3 


N2(2,l) 

— 

6 

N2 (2 , 2) 

= 

8 

AKB(l) 

= 

412.85503 

AKB(2) 

= 

773.599673 

AKS(l) 

= 

6370000.0 

AKS (2) 

= 

773.599673 

T0L1 (1) 

= 

0.01 

TOL1 (2) 

= 

0.01 

TOL6 (1 ) 

= 

0.1 

TOL6 (2) 


0.01 


CK(1 , 1) 


5971137.6 


CK(2 , 1) 
CK(2,2) 
El (1,1) 
El (2,1) 
El (2,2) 


T(I) 

A(I) 

184.99481 

10.787315 

412.5503 

8.92405155 

523.87973 

10" 20 

528.87973 

1.96132994 

800.0 

1.96132994 


XO 

YO 

XDO 

YDO 

SGO 

RO 


6552000.0 


1.5707963 


0.1 


0.1 


io- 5 


B(I) 

C(I) 

1.2878664 

-0.261282644 

0.69140715 

-0.212616813 

1.0 

0.0 

0.36021199 

-0.046948359 

0.36021199 

-0.046948359 

236255.46 


6456667.82 


2800.5746 


635.94317 


9.81 


6370000.0 





TABLE 1 


COMPUTER REQUIREMENTS AND 

Floating point 
32,000 storage locations 
FORTRAN IV compiler 
Subroutine library 

Job Number 

Estimated running time 
SS 2 On 
SS 2 Off 
Restart procedure 
Deck sequence 
Restrictions 

Printed diagnostics 
Originator 


USER INFORMATION 

Hardware 

Hardware 

Software 

Software 

560490 

0.176 hours 

0.075 hour 

No capability 

See computer Listings 

See Introduction and 
Deck Setup Instructions 

None 

Lyle R. Dickey, 

Technical & Scientific Staff, 
Aero-Astrodynamics Laboratory 
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TABLE 2 


PROGRAM INPUT DATA 


Formula 

Notation 

Computer 

Language Notation 

Units of 
Measurement 

Remarks 

N 

N 



Number of staging times 

m 2 

NO 



Number of options 

% 2 

NR 



Number of max. constraints 

q 

NC 



NC#0 if equation B is used 


10 



Number for the weighting 
function option 


NRUN 



Number of cases 

E o 

EO 



Tolerance for E(t^) 


T2 

seconds 


Desired t 2 value 

p 

N6(I) 



Numbers of power cutoff options 

g 

N5(I) 



Numbers of coast cutoff options 

. 

N2 (I , J) 



Numbers of additional constraints 

k 

Pm 

AKB(I) 

see Tbl 

3 

Values of power cutoff limits 

k* 

K Sm 

AKS(I) 

see Tbl 

3 

Values of coast cutoff limits 

t c 

DTA(I) 

seconds 


Integration increment for 
coast 

TL 

TOLL (I) 

see Tbl 

3 

Tolerance for power cutoff 

TL“ 

T0L6 (I) 

see Tbl 

3 

Tolerance for coast cutoff 

K 

CK(I, J) 

see Tbl 

3 

Values of cutoff constraint 
limits 
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TABLE 2 

PROGRAM INPUT DATA 
(Continued) 


Formula 

Notation 

Computer 

Language Notation 

Units of 
Measurement 

Remarks 

EL l 

El (I) 

see Tbl 3 

Tolerance for AZ 

T i 

T(I) 

seconds 

Staging times 

H 

A(I) 

Newtons 

Constant times force 

b i 

B(I) 

mass 

Constant times mass 

c i 

C(I) 

mass/ seconds 

Constant times m 

x o 

XO 

meters 

Horizontal initial distance 

y Q 

YO 

meters 

Vertical initial distance 

x o 

XDO 

meters/second 

Horizontal initial velocity 

y Q 

YDO 

meters/second 

Vertical initial velocity 

§o 

SGO 

2 

meters/ second 

Initial gravity acceleration 

r o 

RO 

meters 

Initial distance from center 
of earth 

At 

DELT 

seconds 

Power integration increment 

% 

THC 

radians 

Velocity vector angle at 
cutoff 


27 





TABLE 4 


DATA PACKAGE CARD ARRANGEMENT 


Card 

Variable 

Data 

Remarks 

Number 


Field 


i 

N 

1-5 

Number of stages + 1 

i 

NO 

6-10 

Number of options 

i 

NR 

11-15 

Number of constraints 

i 

NC 

16-20 

0 if equation B is not used 

i 

NRUN 

21-25 

0 for last case in a stack 

i 

10 

26-30 

Number of weighting function option 

i 

EO 

31-45 

Tolerance for G(tj^) 

i 

T2 

46-60 

Value of t2 desired 

2 

N6(l) 

1-5 


2 

N6(2) 

6-10 

Numbers for power 

• 

• 

• 

Cutoff options, 

• 

• 

♦ 

I = 1,2, 3,---, NO 

2 

N6 (I) 

• 


3 

N5 (1) 

1-5 


3 

N5(2) 

6-10 

Numbers for coast 

• 

• 

« 

Cutoff options, 

* 

• 

• 

-r _ i o o vrr\ 

_L — J- , J 1NW 

3 

N5(I) 

• 
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TABLE 4 


DATA PACKAGE CARD ARRANGEMENT 
(Continued) 


Card 

Number 

Variable 

Data 

Field 

Remarks 

4 

N2(l,l) 

1-5 


4 

N2(l j2) 

6-10 

Numbers for 

• 

• 

• 

additional constraints 

• 

• 

• 

I = 1,2,3,— , NO 

4 

N2 (1 , J) 

• 


5 

N2 (2 , 1) 

1-5 

J - 1,2,3, , NR 

5 

• 

N2 (2 , 2) 

• 

6-10 

• 


• 

NO+3 

• 

N2(I,J) 

• 

• 


NO+4 

AKB(l) 

1-15 

First power cutoff value 

NO+4 

AKS(l) 

16-30 

First coast cutoff value 

NO+4 

DTA(l) 

31-45 

First coast integration increment 

NO+4 

TOL1 (1) 

46-60 

First power cutoff tolerance 

NO+4 

TOL6 (1) 

61-75 

First coast cutoff tolerance 

N0+5 

AKB(2) 

1-15 

Second power cutoff value 

NO+5 

AKS (2) 

16-30 

Second coast cutoff value 

• 

• 

• 

• 

• 

• 

« 

• 
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TABLE 4 


DATA PACKAGE CARD ARRANGEMENT 
(Continued) 


Card 

Number 

Variable 

Data 

Field 

Remarks 

2N0+3 

T0L6 (I) 

61-75 

I = 1,2,3,— , NO 

2N0+4 

CK(1 , 1) 

1-15 

First constraint value for first option 

2N0+4 

El (1,1) 

16-30 

Tolerance for first constraint value 

2N0+5 

GK(1 , 2) 

1-15 

Second constraint value for first option 

2N0+5 

El (1,2) 

16-30 

Tolerance for second constraint value 

2N0+6 

• 

CK(1 , J) 

• 

1-15 

• 

Third constraint value for first option 

• 

• 

2N0+LC* 

• 

El (I, J) 

• 

16-30 

♦ 

Tolerance for last constraint value 

2N0+LC+1 

T (1) 

1-15 

Initial time 

2N0+LC+1 

A(l) 

16-30 

Initial force constant 

2N0+LC+1 

B(l) 

31-45 

Initial mass constant 

2N0+LC+1 

c(D 

46-60 

Initial m constant 


2N0+N+LC 

T (N) 

1-15 

Calculated 

2N0+N+LC 

a 

— \~ ’ / 

16-30 

AfN-1 ^ 

" N" “/ 

2N0+N+LC 

B(N) 

31-45 

B(N-l) 

2N0+N+LC 

C(N) 

46-60 

C(N-l) 



LC = total number 
for one case. 

of constraints 

used 
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TABLE 4 


DATA PACKAGE CARD ARRANGEMENT 


Card 

Number 

Variable 

Data 

Field 

2N0+N+LC+1 

XO 

1-15 

2N0+N+LC+1 

YO 

16-30 

2N0+N+LC+1 

XDO 

31-45 

2N0+N+LC+1 

YDO 

46-60 

2N0+N+LC+2 

SGO 

1-15 

2N0+N+LC+2 

RO 

16-30 

2N0+N+LC+2 

DELT 

31-45 

2N0+N+LC+2 

THC 

46-60 


(Continued) 

Remarks 

Initial horizontal distance 
Initial vertical distance 
Initial horizontal velocity 
Initial vertical velocity 
Gravity at initial time 
Initial radius 
Integration increment 
0 for power cutoff at T(N) 
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TABLE 5 


Notation 

TIME 

DELT 

ETAT 

DETT 

CHI 

DCHI 

X 

Y 

XD 

YD 

XDD 

YDD 

RAD 

VEL 

THETA 

PHI 

OVEL 

EVEL 


OUTPUT NOTATION 
Definition 



A ^k 

^k 

Mt k 

X 

y 

X 

y 

X 

• • 

y 

r 

v 

a 

a> 

°v 


Units 

Seconds 

Seconds 

Seconds 

Seconds 

Degrees 

Degrees 

Meters 

Meters 

Meters/Second 

Me ter s/ Second 

Meter s/ Second 

o 

Meters/Second^ 

Meters 

Meters/Second 

Degrees 

Degrees 

2 

(Meters/Second) 

(Meters /Second) 2 
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TABLE 5 


OUTPUT NOTATION 
(Continued) 


Notation 

Definition 

Units 

ENEG 

E 

Mas^Seconds 

ENG 

E. 

3 

Meters 

HB1 


V Seconds 2 

HB2 

h . 

tC- 

^Seconds 2 

KB1 

k 

l 

X /Seconds 2 

KB2 

K 

CL 

^Seconds 2 

FOM 

f 

m 

Meters/Second 2 

HS1 

h i 

1 /Radians 

HS2 

a 

1 /Radians 2 

H2S 

JL. 

h“ 

2 

X /Radians 2 

EPS 

6 k 

1 /Radians 

DRAGX 

• • 

"g 

Meters /Second 2 

DRAGY 

y'g 

Meters/Second 2 


APPENDIX I 


WRITTEN BY : 
PURPOSE : 

DESCRIPTION : 
USAGE : 


RUNKUT 


Runge-Kutte Integration 
INO 


C. R. Wallace 

To numerically integrate a set of first and second order 
differential equations . 

This FORTRAN subroutine uses the fourth order Runge-Kutta 
integration formulas (Hildebrand, page 238) to numerically 
integrate a set of differential equations which may be 
solved for the highest order derivatives of each variable. 
There must be only one independent variable. 

The subroutine must be called four times for each time 
step. The calling sequence is: 

CALL RUNKUT (NEQ ,NSO,H,X, F, Y, YP,FK, AY, AYP, SUMl , SUM2 , IND) 

NEQ - Number of equations 

NSO - Number of second order equations 

H - Step Size 

X - Independent variable 

F - Highest order derivatives (a block of NEQ values) 
Y - Dependent variables (a block of NEQ values) 

YP - First derivatives of dependent variables (second 
order equations only) (NSO values) 

FK, AY , AYP, SUMl , SUM2 - Auxiliary block storage (NEQ 

values of FK,AY, SUM2; NEQ values 
of AYP and SUMl) 

IND - Indicator for which of the four passes is being 
made through the subroutines; i.e., 

IND = 1,2,3, or 4. 

In order to have all values (including any variables com- 
puted from intermediate equations) correct at each time step 
it is recommended that subroutine RUNKUT be called before 
the highest order derivatives are computed. 
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Example : 


Suppose the equation 


RESTRICTIONS: 


TIMING: 


u" = g L 

(u\ 

U, 

V, 

w 1 , 

w. 

t) 

v' = g 2 

(u‘, 

U, 

V, 

W ! , 

W, 

t) 

w" = g 3 

(u'. 

U, 

V, 

w' , 

w, 

t) 


3 

2 

5 

4 


6 

1 


are to be integrated for 100 steps with a step size 
of .1. This could be coded as follows: 


• 

EQUIVALENCE (u", F(l)), (w", 
(u, Y(l) ) , (w, 

(u', YP (1 ) ) 5 (w 1 , 
DIMENSION F (3) , Y(3), YP(2), 
SUM1 (3) , SUM2 (3) 

NEQ = 3 
NSO = 2 
H = .1 
X = 0. 

NP = 101 
DO 1 I = 1, NP 
IF (1-1) 2, 3, 2 
IND = 4 


F (2) ) , (v' , F (3) ) , 

Y (2) ) , (v, Y (3) ) , 

YP(2) ) 

FK(3) , AY (3) , AYP (3) , 


GO TO 4 
IND = 1 

CALL RUNKUT (NEW ,NS0,H,X,F , Y, YP,FK , AY, AYP, SUMl , SUM2, IND) 


u" = g l 

v ' = g2 

w" = g 3 

IND = IND + 1 

IF (IND - 4) 5, 5, 6 


(any desired coding) 


CONTINUE 


Equations may not be higher than second order. If there 
are both first and second order equations, the variables 
from the second order equations should be placed first in 
the blocks. The error is of the order of times the 
fifth derivative of the dependent variable. 

The computation within the subroutine would require less 
than one second for 4 equations for 100 Lime steps. 
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APPENDIX II 


SUBROUTINE: 
WRITTEN BY: 
MODIFIED BY 
PURPOSE: 
USAGE : 


Matrix Inverse 


Sylvia Bryant 

Compute the inverse of a matrix 
CALL MINV (A, AI, DET, M) 


Where: A is a (MXM) matrix 

AI is the resulting inverse of A 

DET is the value of the determinant of A. If A is 
a singular matrix, DET = 0. 

M is the size of A. 


There exist two versions of this subroutine. In one case A and 
AI are two-dimensional arrays (M, M) and must be dimensioned in the 
user’s program and in the subroutine. In the second version A and AI 
are single arrays and must be dimensioned by the user’s program only. 

In both cases the array I0RD must be dimensioned (M) in the 
subroutine. 


Upon return from the subroutine, A is unchanged. 
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APPENDIX III 


SUBROUTINE: 
WRITTEN BY: 
LANGUAGE: 
PURPOSE: 

USAGE: 

WHERE: 


RESTRICTIONS: 


Matrix Multiply 

James W. Hilliard and Sylvia Bryant 
FORTRAN IV 


To compute the product of (single) (single) , 

(double) (double) , (single) (double) , or 
(double) (single) array matrices. 

CALL MATMP(I1, 12, 13, 14, 15, 16, I, J, K, A, B, C) 


11 is the number of dimensioned 

12 is the number of dimensioned 

13 is the number of dimensioned 

14 is the number of dimensioned 

15 is the number of dimensioned 

16 is the number of dimensioned 
A is a (I x J) matrix 

B is a (J x K) matrix 
C is a (I x K) matrix 


rows of A 
columns of A 
rows of B 
columns of B 
rows of C 
columns of C 


(1) II, 12, 13, 14, 15, and 16 must NOT be zero (0) 
for single array matrices. 


EXAMPLE: DIMENSION A(20), II or 12 must = 1. 

(2) Neither A nor B may be C. 

EXAMPLE: CALL MATMP(I1, 12, 13, 14, 15, 16, 

I, J, K, A, B, A) is illegal. 

(3) I < II J < 12 K < 14 

I < 15 J < 13 K < 16 


STORAGE: 


263 8 - 179 10 


locations . 



FLOW CHART 


FOR COMPUTER PROGRAM 
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NO 


I = 1 


YES 
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I 



r\ n 


m 


DD 50 J = 1, LA 



47 


= 1 
0.0 












CONTINUE 
ON DO 78 
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NUMBERS VALUES TOLERANCES 
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LISTINGS FOR THE COMPUTER 


PROGRAM AND THE THREE SUBPROGRAMS 



DIMENSION CM ( 1 3 ,4 ) >PSM ( 1*4) »CXDM<4*1 ) ,PR(8*4)»CR<8*4).AID<4*4)* 

10M ( 1 *4) *ULM(4,4) »Hl ( 4 ♦ 1 ) »H2 ( 4 * 1 ) * FM ( 8 * 175 ) »GM ( 8 ♦ 175 ) *CAP(8* 1 ) * 

2BM ( 8 » 8 ) >FK (2) ,AY< 2 > ,AYP< 2 ) »SM1 (2 ) ,SM2( 2 ) »AEM4,4) ,AEK1 (4*4) >PRD1 
3(4*4) »PRP2<4*4)*ZZD<8»1) *ULP<4*4) »QBM( 8*1 ) *F1(8*1)*F2(8*1)*CQ(1*4) 
4»FT ( 1»8)*EN(8»1) *CM 13 ♦ 8 ) *E1(13*8) ♦ CKK ( 1 3 ) . E 1 1 ( 1 3 ) *N5< 13) *JPP( 8 ) 
COMMON N2 < 13*8 ) * AKB ( 8 ) *AKS(8)*DTA(8).T(10)*A(10)»B(10)*C(10). 

1ST ( 175) *SDT ( 175) *ET < 175 ) »DET ( 175 ) *CHI ( 175) »D-CHI (175)*R1(175)*V1(17 
25) *PHI ( 175 ) »OV( 175 ) »EV< 175 ) *FOM( 1 75 ) . HB1 < 1 75 ) *HB2 ( 1 75 ) » AK 1 ( 1 75 ) * 
3AK2 (175)* THET(175)*X1(175)*Y1(175)»XD1(175)* YD 1(175) .XDD1 ( 175 ) * 
4YDD1 ( 175) .21 ( 13) *Z2 ( 13) *N6< 13 ) .TOL1 ( 13 ) *T0L6( 13 ) »EP ( 1*175 ) *H2S( 1 *1 
575) *W1 (2) *W2(2)*W3(2) .ENG (175) . FN J ( 1 7 5 ) . I WAS ( 8 ) ♦ I T ( 8 > *HS1 ( 1 75 > . 
6HS2 ( 175) 

COMMON TOL3 ( 1 3 ) *TOL7 ( 13 ) . AKS1 (8)*AKB1(8)*N3(8),N7(8).DTC<8)» 

IN 1 ( 1 3 » 8 ) 

EQUIVALENCE (XDD.W1 ( 1 ) ) ♦ ( YDD.W1 ( 2 > > » < X . W2 < 1 > ) ♦ < Y » W2 < 2 > ) ♦ 

1 ( XD*W3 < 1 ) ) » ( YD»W3 < 2 ) ) 

1 READ ( 60,2 )N»NO,NR»NC* IR. I 0 » NRUN . EO . T2 

2 FORMAT (7I5.2F15.0) 

READ( 60,201 )( N6(J) *J=1 , 13 ) 

'201 FORMAT (1315) 

READ(60,201)(N5(J)*J=1*13) 

DO 300 I = 1 * NO 

300 READ! 60,3 ) (N2 ( I ,J ) *J=1 *8 ) 

3 FORMAT (815) 

DO 4 1=1, NO 

4 READ( 60,5 ) AKB ( I ) ,AKS( I ) ,DTA ( I ) ,T0L1 ( I ) *T0L6 ( I ) 

5 FORMAT ( 5E15.0 ) 

DO 255 1=1,8 
JPP ( I ) =0 

DO 255 J= 1 . 1 3 
255 CK( J, I )=0.0 
DO 253 J= 1 * NO 
K = 1 

DO 253 1=1,13 
I F ( I -N2(J,K) >253,254,253 
254 READ! 60,501 )ADC1 ,ETL1 
CM I » J ) =ADC 1 
El ( I , J)=ETL1 
K = K+ 1 

253 CONTINUE 

501 FORMAT ( 2F1 5.0 ) 

DO 6 1=1, N 

6 READ(60,599)T ( I ) »A( I ) ,B( I ) ,C( I ) 

599 FORMAT ( 4E1 5 • 0 ) 

READ! 60,7 ) XO , YO * XDO , YDO , SGO , RO ,DELT ,THC 

7 FORMAT (4E15.0/4E15.0) 

NN1=N— 1 

WR I TE ( 61,1 1 ) NN1 ,NO 

11 FORMAT ( 1H1 *////46X ,24H*** OPTIMUM GUIDANCE ***//48X * 12HVEH ICLE HAS 
1 12, 7H STAGES/ 50X.4HAND I2.10H OPTION(S)) 

WR I TE ( 61,12 ) 

12 FORMAT ( 11X, 15H K BAR 11X.15H K STAR 11X,15H TL B 

1AR 11X.15H TL STAR ) 

DO 13 1=1, NO 

13 WRITE(61,14)AKB( I ) ,AKS( I ) *T0L1( I ) ,TOL6( I ) 
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14 F0RMAT(4(11X*E15,8) ) 

WRI TE ( 61 » 1 5 ) 

15 FORMAT ( //11X»15H STAGING TIME 11X.15H All) 11X.15H 

IB ( I ) 11X.15H C ( I ) ) 

DO 16 1=1, N 

16 WRITE (61, 14)T { I ) *A( I ) *B( I ) * C < I ) 

<■ WRITE (61*17)XO»YO»X DO >YDO»SGO»RO*THC*DELT 

17 FORMAT! //5X,6HXO= E15.8 ,5X ,6HY0= E 1 5 . 8 , 5X , 6HXD0= E15.8*5X,6HY 

1D0= E15.8/5X ,6HGO= E 1 5 .8 , 5X ,6hRO= El 5 . 8 , 5 X , 6 hTHC= E15.8,5X*6 

2HDELT = El 5 • 8 ) 

WRITE (61, 124) 

. 124 FORMAT ( //16X»10HCONSTRAINT11X ,15H CONSTRAINT 11X,15H CONSTRAIN 

IT, / 1 6X , 1 OH NUMBERS 11X.15H VALUES 11X.12H TOLERANCES) 

DO 125 J= 1 » 8 
NADD=0 

DO 125 1=1,13 
NADD=NADD+1 

IF C CK ( I » J ) ) 126,125,126 

126 WRITE! 61,127 >NADD»CK( I,J),E1(I,J) 

125 CONTINUE 

127 FORMAT! 16X , 1 5 , 5X , 2 ( 1 IX , E 1 5. 8 ) ) 

RA=57, 2957791 
G01=-SG0*R0**2 
TPN=T(N)-(100»0*B(N) )/C(N) 

CON=( XDO*YO-YDO*XO)/< XO*XDO+YO#YDO > 

CH 1 0= AT AN ( CON ) 

I F ( CH 1 0 ) 1 8 , 19 » 19 

• 18 CHI0=CHI0+3. 1415926 

19 THD = THC-CH 1 0 
K = 2 

I FST=4 
DO 20 1=1, N 
I F ( I.EQ.l )21,22 
21 ST ( 1 ) =T ( 1 ) 

DT=T(2)-T(1) 

TLT=T2-T( 1 ) 

IF(DT-TLT >23,24,25 

23 ST ( 5 ) =T2 

SDT ( 4 ) =T2-T ( 2 ) 

GO TO 26 

24 ST ( 5 ) =T ( 2 )+DELT 
SDT ( 4 ) =PELT 

GO TO 26 

25 IF(DT.LE.DELT)27,28 
27 ST ( 5 ) =T2+DELT 

SDT(4)=ST(5)-T(2) 

26 K=K+3 
CHK=0.0 

- ST (2) =ST ( 1 )+0»5*DT 

ST ( 3 ) =T ( 2 ) 

ST ( 4 ) =T ( 2 ) 

SDT ( 1 ) =0. 5*DT 
SDT ( 2 ) =SDT ( 1 ) 

. SDT ( 3 ) =0 • 0 

I FST=0 
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GO TO 137 
28 ST ( 2 ) = T 2 
SDT ( 1 ) =TLT 
GO TO 137 
22 DT=T ( I ) -T ( 1-1 ) 

I F ( IFST.EQ.O) 138,139 

138 I FST=4 
GO TO 20 

139 I F ( I • EQ • 2 ) 140,141 

141 IF(DT.LE.DELT) 142,143 

143 SDT (K)=DELT-CHK 

ST ( K+l ) =ST ( K ) +SDT ( K ) 

GO TO 144 

140 SDT ( K ) =DELT 

ST ( K+l )=ST (K)+DELT 

144 CHK = T ( I )-ST (K + l ) 

I F ( CHK-DELT >145,225,146 

146 K=K+1 

GO TO 140 

145 SDT ( K+ 1 ) =CH K 
226 SDT ( K + 2 ) = 0 » 0 

ST ( K+2 ) =T ( I ) 

ST ( K + 3 ) =T ( I ) 

K = K + 3 
GO TO 137 
225 SDT ( K+ 1 ) =CH K 
CHK = 0 • 0 
GO TO 226 

142 DOT = DT +CHK 

I F ( DOT-DELT >147,148,149 

147 CHK = T ( I ) -ST ( K — 2 ) 

GO TO 150 

148 CHK = 0 *0 
GO TO 150 

149 CHK = T ( I )-(ST ( K — 2 )+DELT ) 

150 K=K+3 

ST ( K-2 ) =ST (K-3)+0»5*DT 
ST ( K-l ) =T ( I ) 

ST ( K ) =T ( I ) 

SDT ( K-3 ) =0 • 5*DT 
SDT ( K-2)=SDT ( K — 3 ) 

SDT ( K-l ) =0.0 
137 K2=K 

K3=K2+1 
20 CONTINUE 

WRITE(61,269)K,K2 
269 FORMAT ( 2 ( 5X , 14) ) 

DO 29 1=1, K2 
I F ( I • EQ . 1 >60,81 
60 ET ( I ) =ST ( I ) 

GO TO 62 

81 ET Cl) =ST <1 — 1 ) +0 • 5*SDT ( I ) 

DET ( I — 1 ) = E T ( I ) -ET ( I - 1 ) 

62 CHI ( I ) = CH I O + THD* ( ( ST ( J ) -ST ( 1 > > / (TPN-ST( 1 ) ) ) 
DCH I ( I ) =0.0 



29 CONTINUE 

DO 202 11=1*4 
DO 202 J J= 1 *4 
IF(I1.EQ.JJ)203*204 

203 A I D ( I I ♦ JJ ) = 1 *0 
GO TO 202 

204 A I D ( I I » JJ ) =0 • 0 
202 CONTINUE 

HI ( 1 » 1 ) =0 «0 
Hl(2»l )=0.0 
H2<1»1)=0.0 
H2 ( 2 • 1 ) = 0 • 0 
AEK ( 1 * 1 ) = 0*0 
AEK( 1 »2 ) =0*0 
AEK ( 1 » 3 ) = 1 • 0 
AEK ( 1 *4 ) =0*0 
AEK(2»1)=0.0 
AEK(2*2)=0.0 
AEK (2*3) =0*0 
AEK ( 2 *4 ) =1 • 0 
AEK<3*3)=0.0 
AEK ( 3 *4 ) =0»0 
AEK(4*3)=0.0 
AEK ( 4 *4 ) =0 »0 
CM ( 1 * 1 ) =0 *0 
CM ( 1 *2 ) =0*0 
CM( 1»3)=0.0 
CM( 1»4 )=j0.0 
DO 75 11=2*5 
DO 75 J J= 1 » 4 

75 CM (II * JJ) = A I D ( I I-1*JJ) 

DO 170 MO=l *NO 
170 I WAS ( MO ) =K2 

30 LI = 1 

I TER= I TER+1 
I DEL=0 
NCOUTT= 1 
NSEC=0 

DO 808 11=1*8 

JPP( II ) =0 

TOL3< I I ) =T0L1 (II) 

TOL7 ( I I ) =T0L6 (II) 

DTC ( I I ) =DTA ( 1 1 ) 

N3( I I ) =N5 ( I I ) 

AKS1 ( I I ) =AKS ( 1 1 ) 

AK8 1 ( I I )*AKB( I I ) 

808 N7( I I )=N6( I I ) 

DO 131 11 = 1* NO 
DO 131 JJ=1 »NR 
131 N 1 ( II *JJ)=N2(II *JJ) 

DO 215 11=1,8 
DO 215 JJ=1 ,175 
FM ( I I * JJ)=0.0 
215 GM ( 1 1 * JJ)=0.0 



LCC = 0 
L2 = 2 
IZ1 = 0 
NC0UNT=0 
I POWER= 1 
I SEC= 1 
LA = NO 
NEQ=2 
NS0=2 

DTR= SDT ( 1 ) 

TRK = ST ( 1 ) 

DO 31 1 = 1 ,K 
I ND= 1 

I F ( I . EQ. 1 ) 3 3 ♦ 34 

33 X=XO 
Y = YO 
XD=XDO 
YD=YDO 
XI ( I ) =XO 
Y 1 ( I ) = YO 
XD1 ( I ) =XDO 
YD1 ( I ) = YDO 

R 1 ( I ) =SQRT(X**2+Y**2 ) 

V 1 ( I )=SQRT(XD**2+YD**2) 

CONl = (X*XD+Y*YD) / ( R1 ( I )*V1 ( I ) ) 

THET ( I ) =ACOS ( CONI ) 

CONl=X/Y 

PHI ( I ) =ATAN ( CONI ) 

OV ( I ) =V1**2+G01/R1 

E V ( I ) =V1**2+2.*G01/R1 

ENG ( I ) =XD*Y-YD*X 

ENJ ( I ) =R1 ( I ) -ENG ( I ) / ( -GOl ) 

FMM=A( I )/B( I ) 

FOM ( I ) =FMM 
CHI I =CHI ( I ) 

XDD1 ( I )=FMM*SIN (CHI I )+X*G01/ ( R 1 ( I )*R1( I) ) 

YDD1 ( I ) =FMM*COS ( CH I I )+Y*G01/ (Rl ( I )*R1 ( I ) > 

34 R=SQRT ( X**2+Y**2 ) 

G=G01/R**2 

DRAGX=X*G/R 

DRAGY=Y*G/R 

XDD=FMM*S I N (CHI I ) +DRAGX 
YDD=FMM*COS(CHI I ) +DRAGY 

CALL RUNKUT (NEQ»NSO»DTR»TRK»Wl ,W2 »W3*FK»AY»AYP »SM1 *SM2»IND) 
IF( IPOWER )38»39*38 
38 IF( IND.EQ. 2)40*41 

40 HB1II ) = (G01/R**3)*(1.-3.*(X/R)**2) 

HB2 ( I ) = { -GOl /R**3)*(3»*X*Y/R**2) 

AKIM ) =HB2 ( I ) 

A K 2 ( I ) = (G01/R**3)*<1.-3.*(Y/R)**2) 

GO TO 47 

41 IF( IND.LT. 3)47*48 

47 AAT = ST ( I )+0.5*SDT( I ) 

I F ( I .EQ.l ) 42 • 24 1 
241 IF(SDT( 1-1 ) .FO. 0.0)42*267 
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267 IF(SDT( I ).EQ. 0.0)46.43 

42 CHI =CH I ( I ) 

CH2=CHI (1+1) 

CH3=CH I ( 1+2 ) 

AAA=AAT-ST( 1+2) 

AAB=AAT-ST( 1+1 ) 

AAC=ST ( I)-ST( 1+2) 

AAD=ST ( I )-ST ( 1 + 1 ) 

AAE = AAT-ST( I ) 

AAF=ST( I +1 ) -ST ( 1+2) 

AAG=ST< 1+1 )-ST( 1 ) 

AAH=ST ( I +2 ) -ST ( I +1 ) 

AA I =ST ( I +2 ) -ST C l ) 

GO TO 44 

43 CH1=CHI ( I -1 ) 

CH2=CHI ( I ) 

CH3=CHI ( 1+1 ) 

AAA= AAT-ST ( 1 + 1 ) 

AAB=AAT-ST ( I ) 

AAC = ST ( 1-1 )-ST( 1 + 1 ) 

AAD=ST ( 1-1 )-ST( I ) 

AAE=AAT-ST ( 1-1 ) 

AAF=ST ( I ) -ST ( 1+1) 

AAG=ST ( I )-ST ( 1-1 ) 

AAH = ST ( 1+1 ) -ST ( I ) 

AA I =ST < I +1 > -ST ( 1-1 ) 

44 CHI I = CH1*< (AAA*AAB)/( AAC*AAD) )+CH2*( (AAA*AAE)/(AAF*AAG) > +CH3* ( (AAB 
1*AAE)/(AAH*AAI ) ) 

45 TAU=( AAT-T(Ll) )/100.0 
FMM=A(L1 )/<B(Ll )+C <L1)*TAU) 

GO TO 46 

48 AAT=ST ( 1+1 ) 

CHI I=CHI <1+1 ) 

CC1=ST ( I )-T(L2) 

IF(CC1.EQ. 0.0)35, 45 

35 IF(ISEC.EQ. 1)37.36 

36 I SEC= 1 
GO TO 45 

37 I SEC=2 
L1=L1+1 
L2=L2+1 
GO TO 45 

46 I ND= I ND+1 

I F ( IND.LE.4>34,49 

49 FOM ( 1+1 ) =FMM 
R=SQRT ( X**2+Y**2 ) 

V=SQRT (XD**2+YD**2 ) 

XDD1 ( I +1 ) =XDD 

YDD1 ( I +1 ) =YDD 
XI ( 1 + 1 )=X 
Y 1 ( I +1 ) =Y 
XD1 ( I + 1 ) =XD 
YD1 ( I +1 ) =YD 
R1 ( 1+1 ) =R 
VI ( 1+1 )=V 
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CON 1 = ( X*XD+Y*YD ) / ( R*V ) 

THET ( 1+1 )=ACOS(CONl ) 

CONl=X/R 

PHI ( 1+1 )=ATAN(CONl ) 

OV( 1+1 >=V**2+G01/R 

EV ( I+l)=v**2+2.*G01/R 

ENG ( 1+1 ) =XD*Y-YD*X 

E N J ( 1+1 )=R-(ENG< 1+1 )**2)/<-G01) 

Z 1 ( 1 ) =ST ( I ) 

zi( 2 )=xim 

Z1 ( 3) =Y1< I) 

Z1 (4) =XD1 ( I ) 

Z 1 ( 5 ) =YD1 (I ) 

Z1(6)=R1U ) 

Z 1 ( 7 ) =V 1 ( I ) 

Z1(8)=THET( I ) 

Z 1 ( 9 ) =PH I (I ) 

Z1 ( 10)=0V(I ) 

Z 1 ( 1 1 ) =EV ( I ) 

Z1 ( 12 )=ENG( I ) 

Z1 ( 13 ) = ENJ< I ) 

GO TO 113 

200 IF(NSEC)113»128*113 
128 Z1 (LY)=Z2 (LY) 

113 Z2( 1 ) = TRK 
Z2 (2 ) =X 
Z2 ( 3 ) = Y 
Z 2 ( 4 ) =XD 
Z2 ( 5) =YD 
Z2 ( 6 ) =R 

Z 2 ( 7 ) = SORT ( XD*XD+YD*YD> 

V=Z 2 ( 7 ) 

CON l=(X*XD+Y*YD)/( R*V ) 

Z 2 ( 8 ) = ACOS ( CONI ) 

CON 1 = X /Y 

Z2(9)=ATAN( CONI ) 

Z2 ( 10 ) =V**2+G01 /R 
Z2( 11 )=V**2+2.*G01/R 
Z2 ( 12 ) =XD*Y-YD*X 
Z2( 13)=ENJ(I+1) 

CM ( 6 » 1 ) =X /R 

CM( 6,2 )=Y/R 

CM( 6,3 ) = 0.0 

CM( 6,4)=0.0 

CM( 7, 1 )=0.0 

CM( 7 , 2 ) =0 • 0 

CM( 7,3 )=XD/R 

CM(7,4)=YD/R 

DMM=R*R*V*SIN(Z2(8) ) 

DMN = R*V*V*SIN(Z2(8) ) 

CM ( 8 , 1 ) = ( X*V*C0S ( Z2 ( 8 ) ) -XD*R ) /DMM 
CM ( 8,2 ) = ( Y*V*C0S(Z2 (8) )-YD*R)/DMM 
CM(8,3)=(XD*R*COS(Z2(8) )-X*V)/DMN 
CMC 8,4 ) = ( YD*R*COS(Z2(8 ) )-Y*V)/DMN 
CM ( 9 , 1 ) =Y /R**2 


76 



CM(9,2)=X/R**2 
CM(9,3)=0.0 
CM(9,4)=0.0 
CM ( 10»l)=GOl*X/R**3 
CM ( 10 * 2 ) =G01 *Y / R**3 
CM( 10 * 3 ) =2 »*XD 
CM ( 10 * 4 ) =2 • *YD 
CM( 1 1 ♦ 1 ) =CM (10*1) 

CM(11»2)=CM(10*2) 

CM ( 1 1 , 3 ) =XD 
CM( 11*4) = YD 
CM( 12 » 1 ) =-YD 
CM( 12 »2 ) =XD 
CM( 12 *3 ) =Y 
CM( 12 »4 ) =-X 

CM(13,l)=X/R+(2.*Z2(12)*YD)/(-G01) 

CM(13*2)=Y/R-(2.*Z2(12)*XD)/<-G01) 

CM( 13*3 )=(-2.*Z2( 12 )*Y)/(-G01 ) 

CM( 13 *4 ) = ( 2 «*Z2 ( 12 ) *X ) / ( -GOl ) 

CXDM( 1*1) =XD 
CXDM< 2*1 )=YD 
CXDM( 3*1) -XDD 
CXDM( 4*1) =YDD 
IF ( I POWER )71*70»71 
71 DO 50 J=1,LA 

IF(JPP(J)“J >256*50, 256 
256 JP=J 

LO=N7 ( JP ) 

TOL2=Z2(LO)-AKBl( JP) 

IF(ABS(TOL2 ) • LE • TOL3 ( JP ) >51*52 
52 DLTC=ABS(Zl(LO)-Z2(LO> ) 

DLTK= ABS ( Z1 ( LO ) -AKB ( JP ) ) 
IF(DLTC.LT.DLTK)83*54 
54 IFCNSEC.EQ. 0)161, 164 
161 I G-K3 

I T ( NCOUTT ) = I +2 
165 SDK IG)=SDT( IG-2) 

ST ( I G ) = ST ( IG-2 ) 

ET ( I G ) =ET ( IG-2 ) 

DET( IG) =DET ( IG-2 ) 

CHI ( IG.)=CHI (IG-2) 

I G= I G-l 

IF < IG.EQ. IT (NCOUTT ) ) 164*165 
164 I DEL= JP 
NStC=4 

151 DZ1=Z2 ( LO )-Zl ( LO ) 

DZ2=AKB1< JP)-Zl(LO) 

CHI ( I+1)=(DZ2*(CHI ( I+D-CHI ( I) > )/DZl+CHI ( I > 
CHI I=CHI ( I ) 

DTR=(DZ2*(TRK-ST< I ) ) )/DZl 
TRK = ST ( I ) 

W2( 1 ) = AY ( 1 ) 

W2 ( 2 ) = AY ( 2 ) 

W3( 1 )=AYP (1 ) 

W3 ( 2 ) = AYR ( 2 ) 


1 = 1-1 
GO TO 31 

83 IF ( NSEC»EQ«.0) 1 88 » 82 
188 DTR=SDT ( 1+1 ) 

GO TO 50 

82 IF ( JP« EG, IDEL ) 151*50 

50 CONTINUE 
GO TO 31 

51 DO 154 M01=l »N 

IF ( TRK.EQ.T (MOl ) ) 155 » 154 

154 CONTINUE 

IFISDTI 1+1 ).EQ.0.0)155*242 
242 IFINCOUTT-l 1227,156*227 

227 IF(AKB(NCOUTT)-AKB( NCOUTT-1 > ) 156 *229,156 
229 I F(N6(NCOUTT)-N6 (NCOUTT-1) ) 156, 155, 156 
156 I WAS ( NCOUTT ) = I 

CHI ( I +2 ) =CH I ( 1+1 ) 

ST ( I + 1 ) =TR£ 

ST ( I +2 ) =TRK 
SDT ( I ) =DTR 
SDT ( 1+1 )=0.0 

SDT ( I+2)=ST ( 1+3) -ST ( 1+2) 

ET ( I +1 ) =ST ( I ) +0 • 5*SDT ( I ) 

ET ( 1+3 ) =ST ( I +1 ) +0 • 5*SDT ( 1+2 ) 

E T ( I +2 ) =ST ( 1 + 1 ) 

DE T ( I ) =ET ( 1 + 1 ) -ET ( I ) 

DET ( I + 1 ) = ET (1+2)— ET(I + 1) 

DET ( I+2)=ET ( I +3 ) — ET ( 1+2) 

DET ( I + 3 ) = FT ( I+4 1-ET ( 1+3 ) 

155 NSEC=0 
JPP( JP)=JP 
I DE L = 0 

GO T0( 58»59)SSWTCHF (2 ) 

58 TH=THET ( 1+1 )*RA 
PH= PH I ( 1 + 1 ) *R A 
CH = CH I I*RA 

WR I TE(61 ,80 ) ITER,TRK,DTR, C T ( I + 1 ) ,DET ( I +1 ) ,CH ,Z2 < 2 ) ,Z2 ( 3 ) ,Z2 ( 4 ) ,Z2 ( 
15)*Z2<6),Z2(7>*TH*PH*Z2(10),Z2<11)»Z2<12>*Z2(13) ,FOM( 1+1 ) »DRAGX ,DR 
2 AGY 

80 FORMAT ( ///32X*39HCONDITIONS FOR POWER CUTOFF, ITERATION I2///5X, 
16HTIME E15.8 ,5X,6HDELT El 5 . 8 , 5X , 6HET AT E 1 5 . 8 , 5X ,6HDETT E15,8/5 
2X ,6HCH I E15#8,5X ,6HX= El 5 . 8 , 5X ,6HY= E 1 5 . 8 , 5X ,6HXD= E15.8 

3/ 5X , 6H YD= E15,8,5X, 6HRAD E 1 5 . 8 , 5X * 6HVEL E 1 5 . 8 , 5X , 6HTHET A E15 

4.8/5X,6HPHI E15.8,5X,6HQVEL E 1 5 . 8 , 5X , 6HEVEL El 5 . 8 , 5X , 6HENEG E 

515»8/5X ,6HENG E 1 5 . 8 , 5X , 6HF0M E 1 5. 8 , 5X , 6HDRAGX E 1 5 . 8 , 5X , 6HDRAGY 

6 E15.8) 

59 T0L2=AKS1 (JP)-AKBl (JP) 

DO 68 11=1,4 

DO 68 J J= 1 * 4 
68 ULM( I I , JJ )=AID( I I ,JJ) 

I F ( N7 ( JP) .EQ.N3 ( JP ) ) 304,86 
304 I F ( ABS( T0L2 ) .LE.0.01 ) 6 1 ,86 
86 I F ( LO-1 )244 ,245 ,244 

245 DO 246 11=1,4 

246 PSM ( 1 ♦ I I )=0.0 
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GO TO 247 
244 DO 238 11=1*4 

238 COt 1. I I ) =CM ( LO* I I ) 

CALL MATMPt 1 » 4 * 4 . 1 * 4 * 4 , 1 ,4 * 1 »CQ * CXDM . PRD1 ) 

DO 239 11=1.4 

239 PSM(1*II)=-CQ(1*II)/PRD1 

247 CALL MATMP<4, 1*1. 4*4*4, 4*1. 4 » CXDM . PSM , PRD1 > 

DO 240 11=1,4 
DO 240 J J = 1 * 4 

240 ULM ( I I » JJ ) = A I D ( I I ,JJ >+PRDl t I I , JJ ) 

GO TO( 53,67)SSWTCHF (4) 

53 WRITE (61, 110) 

110 FORMAT! ///36X,20HCONDITIONS FOR COAST) 

67 I POWER=0 
DTR=DTC ( JP ) 

313 I ND= 1 
FMM=0 • 0 
GO TO 34 

39 I F ( IND.EQ.2 )63,64 

63 HBA=(G01/R**3>*( l.-3.*(X/R)**2 ) 

HBB=(-G01/R**3)*(3.*X*Y/R**2 ) 

HBC=HBB 

HBD=(G01/R**3)*( 1.-3.* ( Y/R)**2> 

64 I ND= I ND+1 

I F ( IND.GT. 4)65*34 

65 AEK (3*1) =HB A 
AEK (3*2) =HBB 
AEK (4*1) =HBC 
AEK ( 4 » 2 ) =HBD 
DO 66 11=1*4 
DO 66 J J= 1 * 4 

66 AEK1( 1 1 , JJ) =AID< 1 1 *JJ)+AEK( I I , JJ)*DTR 

CALL MATMP (4,4*4 *4 *4*4 *4,4*4 * AEK1 *ULM*PRD1) 

DO 69 11=1,4 
DO 69 JJ=1 » 4 

69 ULM< I I , JJ ) =PRD1 ( 1 1 *JJ) 

LY=N3( JP) 

GO TO 200 

70 IFILY.EQ. 1)305*306 

305 T0L2=TRK-AKS1 ( JP) 

I F ( ABS ( TOL2 ) .LE.TOL7( JP ) >91 *311 

311 IF(TOL2)67,67,312 

312 DTR=DTC( JP ) - ( TRK-AKS 1 ( JP ) ) 

TRK=TRK-DTC( JP) 

GO TO 307 

306 T0L2=Z2(LY)-AKS1( JP) 

IF(ABS(TOL2).LE.TOL7( JP) 191,73 

73 DLTC=ABS(Z1 (LY)-Z2 ( LY) ) 

DLT K= ABS ( Z 1 ( L Y ) -AKS i LY ) j 
IF(DLTC.LT.DLTK)89,74 
89 GO TO< 103,111 )SSWTCHF(4) 

103 TH=Z2 ( 8 ) *RA 
PH=Z2 ( 9 ) *RA 

WRITE(61,114)TRK,DTR,Z2(2)*Z2(3)*Z2(4)*Z2(5)*Z2(6)*Z2(7),TH,PH, 
1Z2(10)*Z?(11),Z2(12)»Z2(13) * DRAGX * DRAGY 
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114 FORMAT! / /5X*6HTIME E15 • 8 »5X »6HDELT El 5 . 8 . 5X , 6HX= E15.8»5X»6HY 

1= E15.8/5X ,6HXD= E 15 . 8 , 5X , 6HYD= E 1 5 . 8 » 5X * 6HRAD E15.8,5X»6 

1HVEL E15.8/5X,6HTHETA E 1 5 . 8 , 5X , 6HPH I E 1 5 . 8 , 5X , 6HOVEL E15.8.5X 

2.6HEVEL E15.8/5X.6HENEG E 1 5 . 8 , 5X , 6HENG E 15. 8 * 5X » 6HDRAGX E15.8, 

35X .6HDRAGY E15.8) 

111 IF (NSEC)90»67 ,90 
74 NSEC=4 

90 DZ1=22(LY ) — Z 1 (LY ) 

DZ2=AKS1 (JP ) — Z 1 (LY ) 

TRK=TRK-DTR 

DTR=DZ2*DTR/DZ1 

307 W2( 1 ) = AY ( 1 ) 

W2 ( 2 ) = AY ( 2 ) 

W3 ( 1 ) = AYP ( 1 ) 

W3( 2 ) = AYP ( 2 ) 

GO TO 313 

91 WRITE(61.112)ITER 

112 FORMAT ( ///32X »39HCONDIT IONS FOR COAST CUTOFF, ITERATION 12) 

TH=Z2 <8)*RA 

PH=Z2 ( 9 ) *RA 

WRI TE(61»114)TRK,DTR,Z2(2) ,Z2(3) *Z2(4) ,Z2(5) ,Z2(6) »Z2(7) »TH*PH, 

1Z2 ( 10) »Z2 (11 ) »Z2 ( 12 ) ,Z2( 13 ) , DRAGX » DRAG Y 
NSEC=0 
GO TO 92 
61 NCOUNT = NO 
LW=LO 
GO TO 72 

92 NCOUNT = NCOUNT + 1 
NCOUT T = NCOUTT + 1 
LW=L Y 

72 IF(LW-1)330,331,330 

331 DO 332 11=1,4 

332 PSM ( 1 , I I ) =0.0 
GO TO 333 

330 DO 76 11=1,4 

76 CQ ( 1*11 ) =CM ( LW , I I ) 

CALL MATMP ( 1,4, 4, 1,4,4, 1,4,1 ,CQ*CXDM,PRD2 ) 

DO 77 11=1,4 

77 PSM(1,II)=-CQ(1»II)/PRD2 

333 IF(NC**1)206, 207,206 

206 DO 208 11=1,4 

208 QM (1,11 ) =PSM( 1,11) 

GO TO 205 

207 DO 210 11=1,4 
210 QM< 1 , 1 I ) =CM ( 7 , I I ) 

205 LC=0 

DO 78 11=1, NR 
LB = N1 ( JP, I I ) 

IF(LB)79*78»79 
79 LC = LC + 1 
IZ1=IZ1+1 
ZZD ( IZ1 )=Z2 ( LB ) 

CKK< IZ1 )=CK (LB, JP ) 

Ell ( IZ1 ) = E 1 <LB,JP 
DO 790 M I = 1 ,4 
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790 CR ( LC * M I ) =CM ( LB *M I ) 

78 CONTINUE 

CALL MATMP ( 4*1 *1 »4*4*4*4»1 »4*CXDM*PSM»PRD1 ) 

DO 87 11 = 1, A 
DO 87 M J= 1 , 4 

87 PRD1( 1 1 ,MJ)=AID( 1 1 ,MJ)+PRD1( II ,MJ) 

CALL MATMP(8*4*4»4»8*4*LC»4»4fCR»PRDl*PR) 

KTIM= 1+1 
JJ=I+1 
NCOT=LCC+l 
LCC=NCOT+LC-l 

WRITE(61,257)JPP(JP) , JP » LC ,LCC »NC0T » JJ »KT I M , I Z 1 
257 FORMAT ( //5X*6HJPP= I4,5X,6HJP= I4,5X,6HLC= I4,5X,6HLCC= 

IX *6HNCOT = I4,5X,6HJJ= I 4 * 5X »6HKT I M= I4,5X,6HIZ1= 14) 

DO 94 I 1 = 1 * KT IM 
H 1 ( 3, 1 )=COS (CHI ( JJ) )*FOM ( JJ ) 

H1(4*1)=-SIN(CHI (JJ) ) *FOM ( J J > 

H2 ( 3 * 1 ) =-S I N ( CH I (JJ) ) *FOM ( J J ) *0 • 5 
H2 ( 4 * 1 ) =-COS (CHI (JJ) )*FOM( JJ >*0.5 
IF(I 1-1)216*221,216 
216 AEK(3,U=HB1(JJ) 

AEK ( 3,2) =HB2 ( JJ ) 

AEK ( 4 * 1 ) = AK 1 ( J J ) 

AEK ( 4 , 2 ) = AK2 ( J J ) 

DO 218 11=1,4 
DO 218 J1 =1 ,4 

218 AEK 1 ( 1 1 , J 1 ) =A ID ( I1,J1) + AEK( I1,J1)*SDT(JJ) 

CALL MATMPI 4, 4 , 4 , 4 ,4 ,4 , 4 , 4 , 4 , ULM , AEK1 ,ULP) 

DO 219 11=1,4 

DO 219 J1 =1 ,4 

219 ULM( 1 1 , J1 ) =ULP ( 1 1 , J 1 ) 

221 CALL MAT MP (8»4»4»4,8»4»LC»4»4,PR*ULM*CR) 

CALL MATMP(8»4»4*1»8*1»LC*4*1»CR*H1»F1 ) 

CALL MATMP(8,4,4,1,8,1,LC,4»1,CR,H2,F2) 

LF = 0 

DO 88 I l = NCOT »LCC 
LF=LF+1 

FM( 1 1 ,JJ) =F1 ( LF,1 ) 

GM ( I 1 , JJ ) =F2 ( LF , 1 ) 

88 GM ( 1 1 , J J ) =F2 ( LF » 1 ) 

IFUP.EQ. 10)85,84 

85 CALL MATMP( 1,4, 4, 4,1 ,4, 1,4,4 ,QM ,ULM,PSM) 

CALL MATMP( 1,4,4, 1 ,8, 1,1,4, 1 ,PSM,H1 , FI ) 

CALL MATMP( 1,4, 4, 1,8,1, 1,4,1 ,PSM,H2,F2 ) 

HS1 ( JJ )=F1 ( 1 ,1 ) 

HS2 ( J J ) =F2 ( 1 , 1 ) 

84 JJ=KT I M-I I 
94 CONTINUE 

IF (NCOUNT.EQ.NO >222,223 
223 I P0WER=4 

TRK=ST ( 1+1 ) 

DTR=SDT ( 1+1 ) 

X=X1 ( 1+1 ) 

Y=Y1 ( 1+1 ) 

XD=XD1 ( 1+1 ) 
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YD=YD1 ( 1+1 ) 

31 CONTINUE 
222 K2=KTIM 

WR 1 TE ( 61 » 259 ) 

259 FORMAT! //5X.6HBM MAT) 

DO 93 11=1. LCC 

DO 95 JJ=1.LCC 
BM ( I I >JJ)=0.0 
DO 96 KK=1»K2 

96 BM ( I I ,JJ)=BM( I I,JJ)+FM< I I ,KK ) *FM ( J J » KK ) *DET ( KK ) 

WR I TE ( 61 * 260 ) BM ( I I » J J ) 

260 FORMAT ( 5X *E15. 8 ) 

95 CONTINUE 

93 CONTINUE 

DO 97 11 = 1, LCC 
QBM ( I I ,1 )=0.0 
DO 99 JJ= 1 » K2 

99 QBM (11*1) =QBM ( I I ,1 ) -FM ( I I , JJ )*HS1 ( JJ)*DET ( JJ) 

97 CONTINUE 

CALL MINVIBM.BM.TSST »LCC) 

CALL MATMP(8,8»8,1,8,1,LCC»LCC,1 .BM.QBM.Fl) 

DO 98 11=1, LCC 

98 FT ( 1 , 1 I ) = F 1 111,1) 

CALL MATMP ( 1 , 8 , 8 , 1 7 5 , 1 , 1 75 , 1 , LCC , K2 , FT , FM , EP ) 

CALL MATMP (1,8, 8, 175, 1,175,1 »LCC »K2 .FT ,GM»H2S ) 

DO 186 11=1, K2 

EP< I I ) = EP ( I I ) +HS1 (II) 

186 H2S ( II ) =H2S ( I I )+HS2(I I ) 

DO 100 11 = 1, LCC 

100 ZZD ( I I ) =CKK ( I I ) -ZZD ( I I ) 

DO 889 11=1 , LCC 

IF(A8S(ZZD( II ) ) .LT.F1K II ) ) 889, 102 
889 CONTINUE 
302 DO 101 JJ=1 »K2 

I F ( ABS ( EP ( J J ) J.LT.EO) 101,102 

101 CONTINUE 
WRITE(61,129) ITER 

129 FORMAT ( 1H1»////39X»25HLAST ITERATION — NUMBER I 2//4 IX , 20HALL TOLE 
1RANCES CHECK) 

104 I PR I NT = 0 
GO TO 105 

102 GO TO( 106,107 ) SSWTCHF ( 3) 

106 I PR I N T = 1 

WRI TE ( 61 ,108 ) I TER 

108 FORMAT ( 1H1»////32X»39HC0NDITI0NS FOR POWER FLIGHT, ITERATION 12) 

105 DO 659 I P = 1 , K2 
TH=THET( IP )*RA 
PH=PH I ( I P ) *RA 
CH=CH I ( I P ) *RA 
DH=DCH I ( I P )*RA 

WRITE(61,660)ST(IP),SDT(IP),ET(IP),DET(IP),CH,DH»X1(IP),Y1(IP), 
1XD1 ( IP) ,YD1 ( IP) , XDD1 ( IP ) , YDDl ( IP) ,R1( IP) ,V1< IP) ,TH,PH. OV (IP). 

1EV< I P) ,ENG< IP) , ENJ ( IP) ,HB1 < IP ) ,HB2 (IP) * A < 1 (IP) » AK2 (IP) ,FOM( IP) , 
3HS1( IP) ,HS2( IP) ,H2S( IP) ,EP( IP) 

659 CONTINUE 



660 FORMAT ( //5X*6HTIME E15 . 8 ,5X .6HDELT E 1 5 . 8 , 5X , 6HETAT E15.8,5X,6HD 
1ETT E15.8/5X ,6HCHI E 15.8 » 5X ,6HDCH I E15 . 8 , 5X , 6HX= E15.8,5X*6 

2HY= E15.8/5X *6HXD= E 15 .8 » 5X ,6HYD = E 1 5 . 8 , 5X .6HXDD E15.8,5X 

3 * 6H YDD E15.8/5X,6HRAD E15.8.5X* 6HVEL E 1 5 . 8 , 5X » 6HTHETA E15.8. 

45X,6HPHI E15.8/5X .6H0VEL E 1 5 . 8 . 5X . 6HEVEL E 1 5 . 8 » 5X , 6HENEG £15. 
58,5X,6HENG E15.8/5X,6HHB1 E15.8.5X, 6HHB2 E 1 5 . 8 , 5X ,6HKB 1 El 

65 • 8 » 5X » 6HKB2 E15.8/5X .6HFOM E 1 5 . 8 » 5X »6HHS1 El 5 .8 ♦ 5X * 6HHS2 

7E15.8.5X.6HH2S E 1 5 . 8/ 5X .6HEPS E15.8) 

WR I TE ( 61 »663 ) 

663 FORMAT( //50X,12HETA K MATRIX//) 

DO 664 11=1,4 

664 WR I TE (61,665) ( ULM ( 1 1 » JJ ) » JJ = 1 ,4) 

665 FORMAT (4(5X,E15.8) ) 

I F ( I PRINT )107»109»107 
.1 07 DO 115 11=1, LCC 
DO 116 JJ=1 , LCC 
BM( I I , JJ ) =0.0 
DO 117 KK=1 ,K2 

117 BM( I I , JJ)=BM( 1 1 , JJ )+FM( II »KK ) * ( -FM ( JJ ,KK ) *DET ( KK ) / (2.*H2S< IOC) ) ) 

116 CONTINUE 

115 CONTINUE 

DO 118 11=1, LCC 
EN( I I ,1)=0.0 
DO 119 KK = 1,K.2 

119 EN( I I ,1 )=EN < I I ,1 ) + <FM< I I ,KK ) *EP ( KK ) *DET ( KK ) / ( 2.*H2S ( ) > ) 

118 CONTINUE 

CALL MINV(BM,BM,TSST »LCC) 

DO 120 11=1, LCC 

120 EN( I I ,1) = EN( I I ,1 )+ZZD( I I ,1) 

CALL MATMP(8»8»8*1»8»1»LCC»LCC»1»BM»£N»ZZD) 

DO 121 11=1, LCC 

121 FT ( 1 , 1 I ) =ZZD ( 11,1) 

DO 122 11=1, LCC 

DO 122 KK=1 ,K2 

122 FM( I I ,KK)=-FM( 1 1 » KK ) / ( 2 • *H2S ( KK. ) ) 

CALL MATMP(1, 8, 8,175, 1,175,1, LCC, K2,FT*FM,DCH I ) 

DO 187 II=1,IC2 

187 DCH I ( 1 1 ) =DCH I ( I I )~EP( 1 1 )/(2.*H2S( II)) 

DO 123 KK=1 »K2 

123 CHI (KK)=CHI (KK) +DCH I ( KK ) 

MOVE=2 

DO 159 M02=1,K 
IF(M02.£Q.IWAS( MOVE ) ) 15 7 , 159 

157 I G=M02+ 1 

DO 158 M01=IG,K 
ST ( MOl ) =ST ( MOl+2 ) 

CHI (MOl )=CHI (MOl+2 ) 

SDT ( MOl-1 )=ST ( MOl ) -ST ( MO 1-1 ) 

ET i MOl — 1 ) =ST i MOl—2 ) +0.5'»'SDT ( MO 1-2 ) 

DET ( MOl-1 )=ST (MOl-l)+0.5*SDT (MOl-l)-ET (MO 1-1) 

158 CONTINUE 
MOVE=MOVE+l 

159 CONTINUE 
K4= I 

DO 168 MO=K4,K3 
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168 CHI ( MO+1 ) =3 •*CHI (M0)~3.*CHI (MO-D+CHI < MO-2 ) 

GO TO 30 

109 IF<NRUN)1*252*1 
? 5 ? STOP 
END 

SUBROUTINE MA TMP ( 1 1 * I 2 » I 3 » I 4 * I 5 * I 6 » I » J »K * A ♦ B » C ) 
DIMENSION A ( 1 ) »B( 1 ) »C( 1 ) 

L-l = 0' 

DO 1 L = 1 * K 
DO 2 N= 1 * I 

I C= I 5* I 6- ( I 5* ( 1 6— LI >-N) 

C ( 10 =0.0 
L2 = 0 

DO 3 M*1 *J 

IA=Il*I2-( I 1*( I2-L2 )-N) 

IB=I3*I4-( I 3# ( I 4- LI ) -M ) 

C ( 10=0 IO+A( IA)*B< IB) 

L2=L2+1 
3 CONTINUE 
2 CONTINUE 
L1=L1+1 
1 CONTINUE 
RETURN 
END 

SUBROUTINE Ml NV ( A * A I »DET » M ) 

DIMENSION A( 8*8 ) *AI ( 8*8) * I0RD(8 ) 

N = M 

DO 100 I = 1 * N 
DO 100 J=1*N 
100 A I (I *J)=A( I *J) 

DO 51 I = 1 * N 
51 I ORD ( I ) = I 
DET = 1 . 

DO 1 K = 1 ♦ N 
KP1 =K+1 

TEST=ABS(AI (K*K) ) 

I ND=K 

I F C K— N ) 21 *22*21 
21 DO 23 L=KP1»N 

TEST1=ABS(AI (L*K) ) 

IF(TEST-TFST1)24*23*23 

24 TEST=TEST 1 
I ND=L 

23 CONTINUE 

I F ( I ND-K )25*22*25 

25 DO 26 J=1*N 
TEST1 = AI ( I ND » J ) 

A I ( IND*J)=AI ( K » J ) 

26 At ( K ♦ J ) =T FST 1 
dft=-det 

LL= IORD(K ) 

IORD(K)=IORD( IND) 

IORD( IND)=LL 
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22 CONTINUE 

IF(AI(KtK) )31»32»31 

31 AI (K*K)=W/AI <K»K> 

DO 2 J»1»N 

I F ( J— K ) 3 ♦ 2 » 3 

3 QUOT=AI(K»J)*AT<K»K) 

AI <K»J)=QUOT 

I F ( QUOT ) 42 » 2 *42 
42 DO 4 I = 1 » N 
IF ( I-K ) 5*4*5 

5 TEST* A I (1»J)-AI(I»K) *QUOT 
IF(ABS(TFST)-ABS(AI ( I»J) )*1 . 0E-5 ) 70 .71 *71 

70 A I ( I»J)*0. 

GO TO 4 

71 AI ( J ♦ J ) *TEST 

4 CONTINUE 
2 CONTINUE 

DO, 11 I *1 » N 
I F I I-K) 12*11.12 
12 AI ( I »K ) *-A I ( I *K )*AI ( K *K ) 

11 CONTINUE 

DET»DET/AI <K»K) 

1 CONTINUE 
55 DO 52 ! = 1*N 
I MD“ IORD 1 1 ) 

IF< IND-I >53*52*53 

53 LL=lORD( I ) 

IORD(J )*IORD( IND,) 

I ORD( IND)=LL 

DO 54 K= 1 * N 
TEST=AI(K.IND> 

AI ( K » I ND ) =A I (K»I) 

54 AI ( K * I ) =TEST 
GO TO 55 

52 CONTINUE 
33 RETURN 

32 DET=0. 

RETURN 

END 

SUBROUTINE RUNKUT < NEQ.NSO *H * X »F . Y . YP . FK . AY » AYP »SUM1 »SUM2 » I NDT ) 
DIMENSION F( 100) *Y( 100) *YP(100) »AY ( 100 ) ♦ AYP ( 100 ) *SUM1( 100) 
DIMENSION SUM2( 100) *FKI 100) 

IFUNDT-2) 1*2*5 
1 ALl=H/2. 

AL2=0. 

AL3=» 5 
AL4=1. 

DO 6 1=1 »NEQ 
A Y ( I ) = Y ( I ) 

AYR ( I ) =YP ( I ) 

SUM1 ( I ) =0« 

6 SUM2 ( I ) =0 , 

GO TO 7 ■ 



X=X+H/2 • 

AL2=H/4. 

AL 3 = • 5 
AL4=2. 

GO TO 7 

5 IF (INDT-3) 3*3*4 

3 AL1=H 
AL2=H/2. 

AL3= 1 • 

AL4 = 2 • 

7 DO 8 1=1 *NEQ 

IF (I -NSO) 9*9*10 
9 Y ( I) = A Y ( I )+ALl*AYP( I )+AL2*FK( I ) 
FK< I ) = H*F ( I ) 

YP( I) = AYP ( I )+AL3*FK( I ) 

5UM1 ( I ) =SUM1 ( I)+FK< I ) 

GO TO 11 

10 FK ( I ) =H*F ( I ) 

Y ( I ) = A Y ( I )+A|_3*FK( I ) 

11 SUM2( I )=SUM2( I )+AL4*FK( I ) 

8 CONTINUE 
RETURN 

4 X=X+H/2. 

DO 12 1=1, NEO 
FK ( I ) =H*F ( I ) 

SUM2U > = ( SUM? ( 1 )+FK( I ) ) /6. 

IF (I-NSO) 13*13*14 

13 Y(I )=AY< I)+H*AYP( I)+H* SUMl(I)/6. 
YP ( I ) = AYP ( I ) +SUM2 ( I ) 

GO TO 12 

14 Y ( I ) =A Y ( I ) +SUM2 ( I ) 

12 CONTINUE 
RETURN 
END 
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